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We implement a code to find the gravitational news at future null infinity by using data from a Cauchy code 
as boundary data for a characteristic code. This technique of Cauchy-characteristic Extraction (CCE) allows 
for the unambiguous extraction of gravitational waves from numerical simulations. We first test the technique 
on non-radiative spacetimes: Minkowski spacetime, perturbations of Minkowski spacetime and static black 
hole spacetimes in various gauges. We show the convergence and limitations of the algorithm and illustrate 
its success in cases where other wave extraction methods fail. We further apply our techniques to a standard 
radiative test case for wave extraction: a linearized Teukolsky wave, presenting our results in comparison to the 
Zerilli technique and we argue for the advantages of our method of extraction. 

PACS numbers: 04.25.Dm, 95.30.Sf, 97.60.Lf 

I. INTRODUCTION 

The importance of gravitational waveform templates for gravitational wave detectors implies a need for accurate 3D numerical 
simulations of isolated sources such as binary black hole mergers. These simulations are often done with Cauchy codes based 
on a "3+1" slicing of spacetime. With the slicing conditions most commonly used in numerical simulations, two problems arise: 
artificial boundary conditions must be placed on the computational domain, and information such as the gravitational news 
cannot be extracted at future null infinity X+ . 

To avoid these problems two possible approaches have been suggested. One way is to evolve the entire spacetime. For 
example, a hyperboloidal slicing of the spacetime d 0. IH allows information to propagate to X"*" in finite time. Conformal 
compactifications, such as those suggested by the conformal field equations |4, 5] or employed in ^{ allow X"*" to be located 
at a finite computational coordinate in a regular way. Another example is characteristic evolution based on a Bondi-Sachs line 
element , which gives a natural description of the radiation-region of spacetime extending to X+. Characteristic numerical codes 
have been used to study tail decay |7], critical phenomena L8^ .,9., .10.. .1 L .IZl . singularity structure I.13i il4t il5„ .16.1 and fluid 
collapse to list just a few examples. 

Unfortunately, characteristic methods suffer from the appearance of caustics in the inner strong field region. The problem of 
caustics can be avoided by evolving the strong field region with a standard Cauchy slicing whilst using the characteristic approach 
for the exterior. This technique of Cauchy-characteristic matching (CCE) has proved successful in numerical evolution of the 
spherically symmetric Klein-Gordon-Einstein field equations |20], for 3D non-linear wave equations f2l'] and for the linearized 
harmonic Einstein system |22]. 

The second way of avoiding problems with standard Cauchy codes is a perturbative approach. The standard way of extracting 
gravitational waves from a numerical simulation is based on perturbation theor y, e ither using the quadrupole formula or first 
order gauge invariant formalisms based on the work of Zerilli and Moncrief 11231 l24ll . The vast majority of waveform templates 
are currently based upon these approaches. Numerical codes able to solve a generalization of the Zerilli equation to a three 
dimensional Cartesian coordinate system and to extract the gravitational signal are reported in 125. 26, 27. 281. Extraction of 
gravitational waves based on the Zerilli-Moncrief formalism in fully three-dimensional simulations are presented also in lEollsoll . 

Perturbative methods have been also used to provide boundary conditions at the outer boundary of the Cauchy grid. This 
approach of Cauchy-perturbative matching has been implemented numerically in |31, 32, 33J. All those works are impressive, 
but discrepancies between the results of a perturbative approach and the full non-linear theory cannot be determined without 
solving the fully non-linear problem, although error estimates have been made 

Of these approaches, CCM has many appealing properties. The characteristic description of the exterior, particularly at 
X+, allows for a natural extraction of gravitational information such as the news. Matching to a standard Cauchy code in the 
interior implies that the methods employed in current 3D numerical simulations may immediately be used. As both Cauchy and 
characteristic approaches have been well tested and commonly used, CCM is a natural way of avoiding the potential problems 
of caustics in the characteristic region and artificial boundaries in the Cauchy region. 

In this work we take a step towards the full CCM method by using a Cauchy code to provide boundary data for a charac- 
teristic code which propagates the solution to X+ to extract the waveform. This procedure of Cauchy-characteristic extraction 
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(CCE) allows the computation of gr avitational waves in an unambiguous fashion. CCE was succesfully implemented in the 
quasispherical approximation in ll35ll : here we demonstrate it in the fully non-linear case. 

In this work the outer region extending to X+ is numerically evolved by the Pitt Null Code iMl3ll36l |3^p^p^E^l4lll42ll. 
The link between the Cauchy and the characteristic modules is done by a non-linear 3D CCE algorithm [35, 41, 42, 43,44]. At 
the outer edge of the characteristic grid the Bondi news is computed [37, 39, 45. 461 . We have imported the CCE code into the 
Cactus computational infrastructure i47l l4^ l4^ l5(ill . Within this infrastructure we have been able to use two separate Cauchy 
codes implementing the BSSN and harmonic 3 + 1 formulations of the Einstein equations. This allows us to test the robustness 
of the CCE approach. 

We compare CCE with the Zerilli extraction technique I^3[ l24ll5lll for both non-radiative and radiative spacetimes. The Zerilli 
formalism is based on the paper of Nagar and Rezzolla 1 52'j which reviews and collects the relevant expressions related to the 
Regge-Wheeler and Zerilli equations for the odd and even-parity perturbations of a Schwarzschild spacetime. The conventions 
presented in their review are implemented in the Wave Extract code, included in the Cactus computational toolkit open source 
infrastructure |49]. 

For non-radiative spacetimes, the comparison gives results consistent with reference i3: the CCE algorithm is ©(A^) 
accurate (where A is the computational grid-step), while the accuracy of the perturbative (Zerilli) approach is 0(r^^) where 
r is the radius of the wave extraction sphere. Since the signal is 0{r^^), this implies an 0{r^^) relative error in the Zerilli 
approach. As a result, it is shown in 13411 that CCM is more efficient computationally in the sense that, as the desired error goes 
to zero, the amount of computation required for CCM becomes negligible compared to a pure Cauchy computation. 

We further apply our algorithms to the study of the propagation of a linearized Teukolsky gravitational wave I27.^53l l54ll5^ 
and we compare the CCE and Zerilli wavesignal. For the case of a large extraction radius, we find that both methods give very 
good results and we demonstrate convergence to the analytical waveform of the Teukolsky solution for the CCE news. We show 
that the CCE waveform does not depend upon the extraction radius, which is a major advantage of the CCE method. A small 
extraction radius has no effect on the CCE waveform but it introduces errors in the Zerilli waveform. 

In Sec. In] we outline some general background and notation. In Sec.|lll]we detail the transformation from the Cauchy slice 
to the characteristic slice. In Sec. lIV Al we give a brief summary of the characteristic evolution code. In Sec. lIVBl we describe 
how the news is computed at X+. Finally, in Sec.|V]we give two sets of tests in full 3D numerical relativity to validate the 
accuracy, convergence and robustness of the CCE algorithm. The first set contains four non-radiative tests, with no gravitational 
wave content. The second one is a radiative three-dimensional test. The results show that CCE is valid and accurate in both 
non-radiative and truly radiative situations. 



II. NOTATION, GEOMETRY AND METRICS 

The implementation of CCE described here follows previous descriptions of Cauchy-characteristic extraction and matching in 
the literature. Much of the work has been presented earlier iJa lill l42l l43l 14411 . Here we briefly outline the notation, geometry 
and metrics used. 

The geometry is described by two separate foliations, neither of which cover the entire spacetime. The Cauchy foliation is 
described using a standard 3 + 1 ADM type metric L56ll . 

ds^ = -{a^ - P,(3')dt^ + 2P,dtdx' + -fijdx'dx^. (2.1) 

Many different formulations can be used, given initial and boundary data, to evolve the 3-metric 7^ . In what follows we 
shall either consider the BSSN formulation |57, 581 as implemented in |591 or the generalized harmonic formulation |60] as 
implemented in the Abigel code |61]. In both cases all interaction between the Cauchy and characteristic foliations will be 
performed in terms of the ADM metric, Eq. ( 12. H . 

The Cauchy slice does not extend to asymptotic infinity. Instead an artificial boundary is placed at \x^\ — U. Within this 
artificial boundary a world-tube V is constructed such that its intersection with any t ~ const. Cauchy slice is defined as a 
Cartesian sphere x^ + + — E?, with angular coordinates labeled hy , A — (2, 3). The world-tube F is then used as the 
inner boundary of a characteristic foliation which uses the standard Bondi-Sachs metric I6^l63ll 

ds^ = - {e^'^'^ - r^hABW^U^^ du^ - 2e^^dudr - 2r^hABU^ dudy^ + r^hABdy^dy^. (2.2) 

Here u labels the outgoing null hypersurfaces, y^ — y^ the angular coordinates (the null rays emanating from the world-tube), 
and r is a radial surface area distance. The angular metric Hab obeys the condition 

det{hAB) = dct{qAB), (2.3) 

where qab is the unit sphere metric. This coordinate system consistently covers the world-tube F and the exterior spacetime as 
long as it is free of caustics. 
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The free variables in the Bondi-Sachs metric are then V, [3, and Hab- The physical interpretation of these variables is that 
hjiB contains the 2 radiative degrees of freedom, e^'^ measures the expansion of the nullcone, and V/ r is the counterpart of the 
Newtonian gravitational potential. The 2+1 decomposition of the intrinsic metric on the r = const, world-tube 

-ii.dfdy^ = -e^^-du^ + r^hABidy"^ - U'^du)(dy^ - U^du) (2.4) 
r 

identifies t^Hab as the intrinsic 2-metric of the u foliation, —U^ as the shift vector and e^^V/r as the square of the lapse. 



III. THE CCE ALGORITHM 

The crucial task of the CCE algorithm is to take Cauchy data given in the ADM form Eq. ( 12. U in a neighborhood of the 
world-tube F and to transform it into boundary data for the Bondi-Sachs metric Eq. ( 12. 2> . Then the Bondi code can use the 
hypersurface equations to evolve the appropriate quantities out to X+ so that the gravitational news may be extracted there. 
Much of the present version of the CCE algorithm has been presented in earlier work i3^l4lll4^l43ll44ll . so in this section we 
will give a brief description, highlighting a few new features. 

In section HH the world-tube was defined as a Cartesian sphere + y"^ + = R^, with angular coordinates labeled by 
y"^ , A = {2,3). In addition to the angular coordinates, we setu — t on the world-tube and choose the fourth coordinate, A to be 
the affine parameter along the radial direction, with A|p = 0. The characteristic cones are constructed such that A is the future 
oriented, outgoing null direction normal to the foliation of F. (In order to avoid a singular Jacobian for the Cauchy-characteristic 
coordinate transformation, we require that the world-tube F be timelike.) The affine parameter A is used because the world-tube 
is not constructed to be a surface of constant Bondi r. 

The choice of the angular coordinates is determined following i64ll by the use of two stereographic patches. The use of two 
patches avoids numerical difficulties in taking derivatives near the poles when standard spherical coordinates are used. The use 
of only two patches may not give the most accurate numerical results; various ways of discretizing the 2-sphere on multiple 
patches are discussed in |65]. The two patches are centered around the North and South poles with the stereographic coordinates 
related to the usual spherical coordinates {9, </>) by 



/ l-cosg / l + cos6' ^„,^ 
$North = \ — ^, Csouth = \ -. ^, (3.1) 

1 + cos & V 1 — COS 

and £, = y'^ + iy^ ^ q + ip, where i — \/— 1 . We also introduce the complex vector on the sphere q-^ = (P/2)((5^ + iS^) and 



^3 

its co-vector qa = {2/P){6'^ + iS^), with P = 1 + = 1 + + p^. The orthogonality condition qAq^ = 2 is satisfied by 
construction. The unit sphere metric corresponding to these coordinates is 



1 / N 4 

QAB = 2 [lAQB + qAqB) = 



1 
1 



(3.2) 



The angular subspace in the Bondi code is treated by use of the eth formalism and spin-weighted quantities. (See 13 9116411 for 
details.) 

The CCE algorithm starts by interpolating the Cauchy metric 7^, lapse a and shift /3* and their spatial derivatives onto the 
world-tube. Time derivatives are computed via backwards finite-differencing done along F, e.g., at t = we write 

idtF)[^] = ^ (3F[^] - 4P[^_i] + + O {{At)') . (3.3) 

Knowledge of the Cauchy metric and its 4-derivative is enough to compute the affine metric fj"'^ as a Taylor series expansion 
around the world-tube 

f,'^^ ^fj^^\r + Xfif\r+0{X'). (3.4) 
Next a second null coordinate system (the Bondi-Sachs system) is introduced 

y'' ^{y\y^,y^), where y^=y^, y^ ^ f = u. (3.5) 
The fourth coordinate j/"* = r is a surface area coordinate, defined by 
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The Bondi-Sachs metric on the extraction world-tube can be computed as 

Vi'r-(^^r^) ■ (3.7) 

Note that the metric on the sphere is unchanged by this coordinate transformation, i.e., -q^^ — i)^^ . Therefore one only needs 
to work with the Jacobian components that correspond to derivatives of r. 

In terms of g'^, the two dimensional metric -q^^ can be encoded into the metric functions 

J = ^q^q^'hAB. K = ^^q'^fhAB ■ (3.8) 
The determinant condition Eq. ( 12. 3> translates into 

isT^ = 1 + J J. (3.9) 

With hAB symmetric and of fixed determinant, there are only two degrees of freedom in the angular metric that are encoded into 
the complex function J. 
From Eq. ( 12. 2> we have 

P = -\\og{-ir)='^\og{r^x). (3.10) 

This quantity is a measure of the expansion of the light rays as they propagate outwards. The CCE and the Bondi codes have 
been implemented with the assumption that r a > (or that /3 is real). 
The radial-angular components rf^ can be represented by 

c/=c/V = ^^; (3.11) 



while the radial-radial component 77'"'' is contained in 



V -r 

W = — ^ . (3.12) 



In addition to these quantities, in L40il the auxiliary variables 

V = QJ=^hAB^cq^q''q^\ (3.13) 

k = (5K^ ^hAB^cq^fq^' + ^^K, (3.14) 
B = g/3-/3,^g^, (3.15) 

have been introduced to eliminate the need to explicitly use second angular derivatives in the Bondi evolution code. The required 
boundary data is J, /3, U, drU, W, v, k, and B. (See Sec. lIV Al ) Notice that once the Bondi-Sachs metric is known on the world- 
tube one can only obtain J, f3, U, W. In order to provide the rest of the necessary boundary data we need the radial derivative of 
the Bondi-Sachs metric 

^^^'^ V - dyP ^ ^ dr dXdy^ ^ + dr dyP ' ' ^^^ ^^ 



With ij^J^ akeady known, the only non-trivial parts in Eq.( l3.16> are the Jacobian terms 



These depend on the second derivatives of the Cauchy metric. In order to avoid possible numerical problems caused by interpo- 
lating second derivatives onto the world-tube, we calculate r by taking centered derivatives of r,\ on the world-tube and we 
calculate ?■ au by backwards differencing in time along the world-tube. The remaining term r aa is calculated using the identity 

Aa = (3.18) 
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and the characteristic equation 

P,r^l{j,rJr-{K.r)^) . (3.19) 
o 

(The right hand side of Eq. ( I3.19> can be computed in terms of J ^ — J^\/r\, which, in turn, can be computed in terms of 
already known quantities.) 

Knowledge of the radial derivative of the Bondi metric is important not only for obtaining {drU)^^ but also because the grid- 
structure of the Bondi code is based on the radial coordinate r, and the extraction world-tube will not, in general, coincide with 
any of these radial grid-points. We need to use, therefore, Taylor series expansions to fill the Bondi gridpoints surrounding F 
with the necessary boundary data. We write, e.g., 

/3 = /3|r + A/3,A+0(A2). (3.20) 

Another problem is the need to provide the auxiliary angular variables, _B, k and v on the world-tube. These have been defined 
as the 9-derivatives of Bondi fields in the frame, while taking angular derivatives on the world-tube amounts to computing 3 
derivatives in the frame. The correction term between the two frames is 

= Si^ - ^Br. (3.21) 

The quantities J, 3 J and J a are known from the interpolated Cauchy data, while J aa is not computed in CCE. Thus one can 
compute 9 J but not 9 J.a- As a consequence we cannot use a simple Taylor series expansion to place v on the Bondi grid to 
0{\^) accuracy. The solution to this problem is to first place J on Bondi grid-points surrounding the world-tube, then compute 
9 J on those points (i.e., compute angular derivatives in the Bondi-Sachs frame), and then calculate 9a9 J on the world-tube by 
use of the neighboring Bondi grid values of 5 J. With v ^ and its radial derivative known on the world-tube, we can then 
use the standard 0{\^) expansion to provide v at Bondi gridpoints surrounding the world-tube. This way we make maximal 
use of the Cauchy data as interpolated onto the world-tube and minimal use of the finite difference 9 algorithm that has its own 
discontinuous ©(A^) error at the edges of the stereographic patches. A similar approach is used for k = 9/v and B = 9/3. 

IV. THE BONDI CODE AND THE NEWS ALGORITHM 
A. The Bondi code 

The inner workings of the Bondi code had been described in detail elsewhere (see i20l l35l l36l l37l Issl l39l l40l l4ll l42il ') so 
here we give only a brief overview of the algorithm. As already stated in Sec. Ull] the variables of the code are J, /?, U, W as 
well as Vjk, B. Out of these the equation for J is the only one to contain a time derivative. For this reason J is updated via an 
evolution stencil that involves two time-levels. The rest of the variables are integrated radially from the world-tube out to X+. 
The integration constants are set by CCE. All of these radial integration equations are first differential order in r except for the 
U equation which contains U^rr- For this reason integrating U requires two constants, which explains the need to provide, at the 
world-tube, U as well as U^r- 

B. The News algorithm 

The calculation of the Bondi news function on X+ is based on the algorithm developed in ll39il with the modifications intro- 
duced in 1 45, 46] (an alternative calculation for the news was recently introduced in |6^). Here we present an overview of the 
algorithm. 

The Bondi-Sachs metric Eq. (|^) in a neighborhood of X+ in {u,x^^l = l/r) coordinates (after multiplying by a conformal 
factor P) has the form 

l^ds^ = 0{f)du^ + 2e^^dudl - 2hABU^dudx'^ + hAndx"^ dx^ , (4.1) 

where the metric variables /3, U^, and Hab have the asymptotic expansions (3 — H + 0{P), — + 0{l), and Hab = 
Hab + I Cab + 0(p). We can always find coordinates (ubjQBjPBjIb) (hereafter referred to as 'inertial' coordinates), and an 

associated conformal metric dsB^ = uPds^ {uj > 0), such that (i) (^g§-^^ is null and affine and points along the null generators 

of I+, (ii) Ib — ujI + 0{P), (iii) the conformal metric in the subspace {ub = const., Is = const.) is the unit sphere metric on 
I+. We fix a null-tetrad on X+ by choosing n° to be affine and to point along the null generators of 1+ and by choosing m^"™''' 
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to be the unit-sphere co variant metric in the 2-diinensional angular subspace. In inertial coordinates {uB,qB,PB,lB) on X+ the 
tetrad has the form 

h" = (1,0,0,0) (4.2) 
r = (0,0,0,1) (4.3) 
fh- = (0,P/2,zP/2,0), (4.4) 

where the tilde denotes a quantity defined with respect to inertial observers. Note that rh^'^m^^ — q^^ as required. We define 
a complex vector i^", analogous to m" (i.e. F'^^F^^ — H^^), adapted to the coordinates used in the characteristic evolution 
code via 



F° = (0,i^^,0), (4.5) 



where 



where is the dyad defined in Sec.|lll| Jo = q^q^ Hab/'2', and Kq = q^q^ Hab /2. F° and to" are related by 

m'' ^e"^u}-^F'' + "ffi''. (4.7) 

The 771° term will not enter the news calculation. 

To calculate the Bondi news one needs to evolve two scalar quantities 5, the phase factor in Eq. (14. 7> . and the conformal 
factor Lo, as well the relations f (fs) between the angular coordinates used in the characteristic evolution and inertial angular 
coordinates, and ub{u, ^b) between the inertial time slicing and the time slicing of the characteristic evolution code. Then 6, 
and ub{u,^b) are evolved using the following ODE's along the null generators ofX+ 

^ = -e-, (4.10) 

where Uq ~ qA L^- Also lo is evolved using the PDE 

9„ log w = -5R ^t7 Slog w + igc7^ (4.11) 

The Bondi news function up to a phase factor of e^^"^ is given by 

N = icj-2g-2H^A^B (^^Q^ ^ _ }_c^^j^^lC ^ 2u;DA[u^-^DB{uje'")]^ , (4.12) 

where Da is the covariant derivative with respect to Hab- The Bondi news is calculated in three steps. First Eq. ( I4.12> is 
evaluated, ignoring the e^^*'' phase factor, as a function of the evolution coordinates (u,^). Then using the relation £,{S,b), 
the news is interpolated onto a fixed inertial angular grid (i.e. N{u,£_b)) and multiplied by the phase factor e"^"* (which is 
only known on the inertial grid). Finally, in post-processing, the news is interpolated in time onto fixed inertial time slices (i.e. 
^{ub,£,b))- Once the news is obtained in inertial coordinates it can be decomposed into spin-weighted spherical harmonics. 



Kq + 1 , _4 / 1 



V. TESTS 



We apply the algorithm described above for two sets of tests: non-radiative spacetimes and radiative spacetimes. In the first, 
non-radiative set of tests (section IV AL we analyze Minkowsy and Schwazchild spacetimes. The Minkowski (Sec. IV ATI or 
small perturbation of Minkowski (Sec. IV A2\ tests are used to show the stability of the code and the errors due to transforming 
between the different coordinate systems and sets of variables. 
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The Schwazchild tests describe a static spherically symmetric black hole in a centered frame fSec IV A"3t and in an oscillating 
frame (Sec. IV A"?l . These tests indicate the accuracy of the code in non-trivial spacetimes. 

The last test (Sec. IV Bt is a standard radiative test: namely a linearized Teukolsky wave, and indicates that CCE is a valid and 
accurate method in truly radiative situations. 

In the non-radiative cases, the extracted gravitational wave signal should vanish identically. As this is true uniformly for all 
points at X+ we simplify the computation for the norms of the gravitational wave signal. When the results of the CCE are shown 
the norm used is the simple L2 norm over all grid points at X+, without any weighting by the area element. In the radiative case, 
we plot the real part of the extracted news waveform (9f versus time. 

The Zerilli Moncrief formalism (see |23. 24. 5 1 1 for the implementation used here) approximates the spacetime as a linearized 
perturbation of a spherically symmetric spacetime. Such perturbations are decomposed in terms of spherical harmonics which are 
then expressed in terms of some basis set. An example using the Regge-Wheeler set is where the even parity metric perturbations 
are decomposed as 



^2 iKY'^'-icd + GVdVcF'™) 



Once the basis functions such as G, h!f \ H2 and b are found from the decomposition, the even parity master function 
can be computed from 



{1-2)1 A [r(A - 2) + 6M] 



(5.2) 



Equivalent definitions in terms of other basis functions, alternative conventions and expressions for the odd parity sector can be 
found in llsl . 

Using this notation where gives the odd parity master function, the asymptotic value for the wave forms is 

h+-ih..=^ J2 i J\ QLit')dt') -2Y'"\e, ^)+0 (^1) . (5.3) 

For the non-radiative cases, when we compute the wavesignal from the Zerilli approximation, we use the norm 

mUf = j {h+ - ih^){h+ - ih^)Un (5.4) 

where d£L is the element of solid angle. The last equality Eq. ( 15. 5> follows from the orthogonality of the spin-weighted tensor 
spherical harmonics. We expect hj^ — ih^ to vanish for all angles so this norm should be identically zero. In what follows we 
shall look at the spherical harmonics up to ^ = to = 7 although in practice there is little contribution from the higher modes. 

For the Teukolsky wave used here we expect all terms except for the Q20 term in Eq. j5.3> to be negligible asymptotically, 
which leads to a simple but tedious computation of the wavesignal. 



A. Non-radiative Spacetime Tests 

1. Minkowski Flat Spacetime 

Minkowski space in standard coordinates is evolved in the harmonic Abigel code, with the metric given analytically on all 
Cauchy slices. The domain extends between G [—10, 10], the extraction world-tube is placed at r = 7, and the simulation is 
run until t = 10. The coarsest simulation used 50^ points for the Cauchy grid and 35^ x 31 points for the characteristic grid. 
Simulations to test convergence scaled all grids by factors of 2. For this test, extraction using the Zerilli method gives results 
that are identically zero, as expected. 

This test indicates the level of numerical round off error in transforming from the Cauchy variables to the Bondi variables on 
the world-tube. As shown in Fig. [J the extracted news is small in all cases but the level of truncation error increases linearly 
with the number of points on the extraction world-tube. Since the construction of the boundary data for the CCE code involves 
derivatives of the Cauchy metric, this sensitivity to noise in the Cauchy data is expected. However, given the extremely low 
amplitude of the noise, it is unlikely that this will cause problems in practical simulations. 
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FIG. 1: The L2 norm of the CCE news for Minkowski space in standard coordinates. This indicates the level of numerical noise introduced by 
transforming variables on the world-tube. In this case the noise is linearly dependent on the number of points on the extraction world-tube. It 
is, however, extremely small. 



2. Random Perturbations around Minkowski Flat Spacetime 

This test is in the spirit of the robust stabihty tests of E2l l67ll68il . The stable evolution and extraction of white noise initial 
data with no frequency dependent growth of the wavesignal is a good indication that the combined evolution and extraction 
codes are stable. The domain extends between x' e [—10, 10], the extraction world-tube is placed at r = 7, and the simulation 
is run until t = 100. The coarsest Cauchy grid has 51^ points and the coarsest characteristic grid 35^ x 31 points. 

The results in Fig. |2] show that the error is independent of the resolution of both Cauchy and characteristic grids. This is 
strong evidence that the CCE code is stable against small perturbations that in a practical run would be induced by numerical 
error. The wavesignal from the Zerilli extraction is about an order of magnitude smaller. This is to be expected as the Zerilli 
extraction approach uses the field values. The error in the CCE news has a jump at the beginning. This arises because on the 
initial characteristic data is set to zero so that it takes some time for the noise to build up on the outgoing null cone. 



3. Schwarzchild Black Hole in a "centered" frame 

To test a simple black hole spacetime we use a Schwarzschild black hole in ingoing Eddington-Finklestein coordinates {i, i*), 
with the line element 

ds' = - (^1 - dP + (^^) didf +1^1 + ^^ df' + dn\ (5.6) 

This is manifestly static in these coordinates. The spacetime is evolved using the BSSN code and excision methods described 
in |69|. The coarsest Cauchy grid has 29^^ points with Ax' = 0AA4, whilst the characteristic grid has 35^ x 31 points. The 
world-tube is at r = 7M. In the Cauchy evolution domain octant symmetry is used. 

The evolution is only performed for a short time (to t — lOOM). Over this timescale we see second order convergence for 
the CCE news until t k, 7M and second order convergence in the waveform from Zerilli extraction until t w 20M, as seen in 
Fig.|3] By varying the location of the world-tube or the Zerilli extraction sphere we can see that the errors come from a variety 
of locations. 

The difference between the times at which the two extraction methods lose convergence is probably due to the greater dif- 
ferencing error seen in the Zerilli extraction. At early times (where finite differencing error dominates) the absolute value of 
the error is larger for Zerilli extraction by a factor w 4, as seen in Fig.|3] At late times (where outer Cauchy boundary errors 
dominate) the error for Zerilli extraction is of the same order as for the CCE extraction. 

The non-convergent errors are probably caused by the boundary conditions on the Cauchy grid which do not satisfy the 
constraints. As in Ii69il we are simply applying Sommerfeld type boundary conditions to all fields. This condition does not a 
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FIG. 2: The L2 norm of the CCE news for random perturbations of Minkowski where the Cauchy slice is evolved using the Abigel code. 
The news displays very slow growth that is independent of the frequency of the initial data. This is a good indication that the CCE method is 
stable. The right hand figure shows that the norm of the waveform extracted from the Cauchy grid using the Zerilli method is about one order 
of magnitude smaller. 



priori satisfy the constraints and is not known to be well-posed. Thus we might expect errors to be induced by the use of these 
boundary conditions. It is likely that constraint satisfying boundary conditions or Cauchy-characteristic matching would solve or 
at least greatly reduce this problem. Due to the complicated pattern of the reflected errors arising from a boundary condition on a 
cubical surface, we have been unable to determine exactly how this error depends on the location of the outer Cauchy boundary. 
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FIG. 3: The scaled L2 norm of the CCE news for a static spherically symmetric black hole where the Cauchy slice is evolved using the BSSN 
formalism. At early times, when the world-tube is causally disconnected from the outer boundary of the Cauchy slice, the CCE news converges 
to zero at second order as it should. Later, when the outer boundary is causally connected to the world-tube, deviations from second order 
convergence appear, indicating the effect of the Cauchy boundary conditions on the extracted wavesignal. The left panel shows the CCE news 
extracted at , whilst the right panel shows the "norm" of the wavesignal extracted by the Zerilli method at r = 7. 
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4. Schwarzchild Black Hole in an oscillating frame 

We use the line element for the standard Schwarzschild black hole in ingoing Eddington-Finklestein coordinates, given in 
Eq. ( 15. 6> . The moving coordinate frame (i, x^) is given by 



(5.7) 



where B"^ are parameters specifying the velocity and b{t) is a simple periodic function turned on after some time t > to', here 
we use 



Kt) 



0, t<to 
sin(tj(i — to))^, t > to. 



(5.8) 



For the test shown here we set — 0.2, B^ — 0.5, B^ — 0.3 for the velocity and uj ~ 2n x 0.05, to — 0.5 for the function b. 

The metric is given analytically on the Cauchy grid to avoid any boundary effects. The coarsest Cauchy grid has 51^ points 
and covers the domain x' G [— lOAf, lOA/]. The coarsest characteristic grid has 35^ x 31 points. The extraction world-tube is 
at r = 7M. The simulation is evolved until t = lOOAf . 

The CCE news converges to zero as it should. However, Zerilli extraction does not give a signal that converges to zero as grid 
resolution is refined. Instead, an erroneous non-trivial signal is computed. 
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FIG. 4: The scaled L2 norm of the CCE news for a static spherically symmetric black hole in an oscillating frame. The CCE news converges 
to zero at second order as it should. 



B. Teukolsky Wave Test 

1. Theoretical part 



The Teukolsky solution ll27ll53LI55r to the linearized Einstein equation is a weak gravitational wave propagating through space, 
and represents one of the most valuable standard test cases for numerical codes that implement wave extraction techniques. 
The general form of the spacetime metric is 



-dt^ + (1 + Afrr)dr^ + 2Bfrerdrd0 



+ 2Bfr4,r sin Odrdc/) + {1 + Cf, 



(1) 



A/efydO' 



+ 2{A-2C)fe^r^sin< 



+ (l + C/l'i+^/rV'sin 



(5.9) 
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FIG. 5: The wavesignal computed using Zerilli extraction for the static spherically symmetric black hole in an oscillating frame. The "norm" 
is computed as described in the text with the radial dependence removed. The left panel shows that as the resolution of both the Cartesian 
grid and extraction sphere are increased the signal does not converge to zero. The right panel shows how this erroneous non-trivial wavesignal 
decays as 0{r~^) as the extraction radius is increased. In the main plot the standard scaling of the signal is used. In the inset the curves are 
scaled by the extraction radius so as to overlay each other for the expected decay rate. 



where the angular functions , corresponding to an Z = 2, m — 0, spin — weight — 2 spherical harmonic, are 

frr — 2 — 3 sin^ 6, fre = —3 sm0cos9, fr4, = 0, 

fo^ = 0, 4y = -4^ - 3 sin^ e-1. (5.10) 

The functions A, B and C are given in terms of a free generating function F{x) — Axe~^^ / X^, where A is the amplitude 
and A determines the width of the wave, by 
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dlF 


3d^F 


--) 




4 

r4 

3dlF 


J.5 ! 

6d^F 


( ;2 + 

r4F ^ 


2dlF 


f 

9dlF 









6F, 



^ ^ + rr^ + rr^ + ^ + (5.11) 



where d"F denotes 



d'^F(x) 

d'^F = , ^ ' . (5.12) 

Here, x — t — r corresponds to an outgoing wave. To obtain the ingoing wave solution coresponding to x = t + r, we change 
the sign in front of all the terms with odd numbers of F. 

We consider a superposition of ingoing t + r and outgoing t — r waves centered at the origin of the coordinate system 
at i = 0, which provides a moment of time symmetry. This solution, after is linearized in amplitude, is implemented in the 
harmonic Abigel code. 

We further give the analytical form of the Teukolsky wave signal. The real part of the Bondi news function is proportional 
to the time derivative of the "plus" polarization mode of the gravitational wave signal. The time derivative of the CCE news 
satisfies 



N = lim r'1'4 



(5.13) 
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where ^4 is the Newman-Penrose component of the Weyl tensor 

^"4 = -Cf,p^rn''fhPn''fh^ (5.14) 

with n"^ an ingoing null vector, and mP a complex unit vector, oriented in the angular directions. 
The connection between and the dyad vector q-^ on the unit sphere is 

= -q"^. (5.15) 

r 

Moreover, n'^df^ = dx- A straightforward computation starting from Eq. ( I5.13> . gives 

N = ^-q^q''dlgAB- (5.16) 
r 

At infinity, in the linearized regime, only the 0(l/r) terms in the Teukolsky functions A, B and C contribute to the news. The 
term that determines the gravitational wave signal is d'^.F. We find 

N = -^-^dtF. (5.17) 
Alternatingly, computation of the hnearized expression for the Bondi news from Eq. (I4.12> . gives 

N = -^-^dtF. (5.18) 



C. Numerical part 

We give the metric specified by the Teukolsky solution analytically on the Cauchy grid. Then data is extracted at the world- 
tube and evolved by the CCE code. We carry out a series of simulations, varying both the location of the characteristic world-tube 
and the radius of the ZerilH sphere. We use SO'^, 120^, and 160"^ gridpoints for the Cartesian Cauchy grid and 60^ x 80, 90^ x 120, 
and 120^ x 160 gridpoints for the null grid. The domain extends between G [—15,15] and the simulations are run until t = 30. 

We study the dependence of the signal with the amplitude and conclude that CCE code resolves correctly amplitudes of 
A > 10^^. At smaller amplitudes, clean convergence behavior is contaminated by round-off error. We show results only for 
A = 10^^, A = 1. We study also the dependence of the wave signal with the world- tube radius and conclude that the accuracy 
of the computed CCE news is preserved even for radii as small as r = 5. 

Fig-Elshows the convergence of the CCE news to the analytical solution, for a world-tube radius of r = 10. The convergence 
rate of the CCE news to the analytical data, is given by 

, / 11^80 - NanaW n, ,c in\ 

Cr=i0g2{j— r), (5.19) 

||iVi60 — JVaiiall 

The convergence rate of the computed CCE news to the analytic value Eq. ( I5.18> at i — r = 0, corresponding to the peak of the 
radiated signal, is 

Cr = 2.159. (5.20) 

Hence, the CCE wavesignal is second order convergent and independent of the world-tube radius, as expected. 

Fig.0shows that for small extraction radii, the Zerilli waveform has a slight asymmetry, which is caused by the dependence of 
the Zerilli formalism upon the extraction radius. We have to increase the extraction radius in order to decrease this error Because 
the error does not fall off sufficiently fast with radius in the near zone, where we can realistically carry out the simulation, we 
instead analytically compute the Zerilli news at the extraction radius r ~ 300. Fig.|8ldemonstrates that the agreement between 
the computed CCE news at small radii and the analytical Zerilli news at big radii, is very good. 



VI. CONCLUSIONS 



We have demonstrated the accuracy of Cauchy-characteristic extraction in 3D numerical relativity. The interface at the ex- 
traction world-tube does not introduce significant error with either the BSSN 1591 or harmonic 123. 16111 formulations. 
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Right panel: the absolute value of the error in Zerilli news with radius (grid resolution of 120'^/90^ x 120. 



The CCE news is stable against small perturbations on the Cauchy grid, as shown in Sec. I V A 21 and the level of truncation 
error caused by transforming from the Cauchy to the characteristic variables is small, as shown in Sec. rVAll In a non-trivial 
black hole spacetime, CCE performs as expected. Second order convergence is found in the moving Schwarzhild black hole 
test in Sec. IV A4l and in the early time behavior of the BSSN test in Sec. IV A 3l However, as seen in the late time behavior in 
Sec. IV A 3l the effect of improper BSSN outer boundary conditions on the Cauchy grid is clearly visible in the extracted CCE 
news. 

In the linearized Teukolsky wave test in Sec. IV Bl we demonstrate that the computed CCE news converges to the analytic 
Teukolsky waveform and does not depend on the world-tube radius. The accuracy of the CCE news is preserved even for small 
radii, while the Zerilli news is affected by near zone error. When the extraction radius is sufficiently large, both Zerilli and CCE 
give excellent results. The advantage of CCE over the Zerilli method is that the extracted CCE news does not depend on the 
extraction radius. 
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FIG. 8; When the extraction radius is sufficiently large, the comparison between the Zerilli and the CCE news is in very good agreement. 

Cauchy-characteristic extraction can be seen either as a first step towards a full Cauchy-characteristic matching code or as 
an improved gravitational wave extraction method in its own right. The results of this paper show that as a stand alone wave 
extraction method, it produces the correct results in situations where other methods, such as Zerilli extraction, may fail. Such a 
situation is seen in Sec. IV A 4l where the ZerilU method fails to converge as a result of the pure gauge motion of the Schwarzchild 
metric. 

Finally, these tests also show the need for improvement of the current implementation of the code. Any angular dependence 
in the solution whether through genuine physics, gauge effects, or due to the imposition of boundary conditions on the cubical 
Cauchy boundary, leads to short wavelenght error that is poorly resolved by the CCE code. This is particularly noticeable for 
features in the region where the stereographic patches overlap. Improvements in this area of the implementation will enhance 
the accuracy of the code in extracting the gravitational news from astrophysical simulations requiring high resolution. To this 
end, we are investigating the use of different multiple patch implementations such as l,65il . 
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